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Abstract 



Multiresolution analyses based upon interpolets, interpolating scaling functions in- 
troduced by Deslauriers and Dubuc, are particularly well-suited to physical applications 
because they allow exact recovery of the multiresolution representation of a function 
from its sample values on a finite set of points in space. We present a detailed study 
of the application of wavelet concepts to physical problems expressed in such bases. 
The manuscript describes algorithms for the associated transforms which, for properly 
constructed grids of variable resolution, compute correctly without having to introduce 
extra grid points. We demonstrate that for the application of local homogeneous oper- 
ators in such bases, the non-standard multiply of Beylkin, Coifman and Rokhlin|Q] also 
proceeds exactly for inhomogeneous grids of appropriate form. To obtain less stringent 
conditions on the grids, we generalize the non-standard multiply so that communica- 
tion may proceed between non-adjacent levels. The manuscript concludes with timing 
comparisons against nai've algorithms and an illustration of the scale-independence of 
the convergence rate of the conjugate gradient solution of Poisson's equation using a 
simple preconditioning. 



1 Introduction and Motivation 

Wavelets offer a means of approximating functions that allows selective refinement. If regions 
of an image or a signal have exceptionally large variations, one need only store a set of 
coefficients, determined by function values in neighborhoods of those regions, in order to 
reproduce these variations accurately. In this way, one can have approximations of functions 
in terms of a basis that has spatially varying resolution. This approach reduces the memory 
storage required to represent functions and may be used for data compression. 

Physical applications often involve multiple physical fields which interact in space with 
non-linear couplings. Efficient implementations must minimize not only storage to represent 
the these fields but also the processing required to describe their interactions. It is highly 
desirable to perform the needed operations with a fixed, limited number of floating point 
operations for each expansion coefficient and to minimize the number of points in space at 
which physical interactions must be evaluated. 

As a concrete example of a realistic application, consider the computation of the quantum 
mechanical electronic structure of a collection of atoms in three dimensions. For other 
examples of physical applications, the reader may wish to consult 1^, et ai. Arias 



and coworkers 0,0, and other works which have appeared after the original submission of 
this manuscript nearly one year ago [^,^6], have studied the use of multiresolution bases in 
quantum mechanical computations. (For a review, see [Q.) It is a consequence of quantum 
physics that, near atomic nuclei, electronic wave functions vary rapidly and, that far from 
the nuclei, the wave functions tend to be much more smooth. From this observation, one can 
anticipate that the fine scale coefficients in a multiresolution analysis of the electronic wave 
functions and associated physical fields will be significant only near the atomic cores, allowing 
for truncation. This problem is thus an ideal candidate for multiresolution techniques. 

Within density functional theory |T^, the quantum physics of electrons and nuclei in- 
volves two types of fields, the Schrodinger wave function {V'j(r)} for each electron i and the 
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electrostatic potential 0(r) arising from the average electronic density n{r) = Z]j|'^j(^)P- 
Within the local density approximation (LDA) pO|, the solution for the correct values of 
these fields is obtained at the saddle-point of lowest energy of the Lagrangian functional 

CLDAimA) = ^E//^'^l|V^iWII' + Jd'rVnnc{rHr) + j dh e,,{n{r))n{r) 

•i 

- I d^r(p{r)n{r) -—I rfV ||V0(r)||2. (1) 

Here, we work in units ^ = m = e = 1, Vnuc('") is the total potential which each electron 
feels due to the presence of the atomic nuclei, and exc{n) is a non-algebraic function known 
only through tabulated values. (For a review of density functional theory, see p3| .) 
In practice, one finds the fields {tpiir)}, (j){r) by 

• expanding the fields in terms of unknown coefficients within some basis set 

a 
a 

• evaluating Eq. (|I]) in terms of the unknown coefficients c and d; 

• determining analytically the gradients of the resulting Clda{c, d) with respect to those 
coefficients; and 

• proceeding with conjugate gradients to locate the saddle point. 

All follows directly once one has expressed the Lagrangian as a function of the expansion 
coefficients. 

In doing this, we note that each term represents a local coupling in space, but that one 
coupling, {/)(r)n(r), is cubic in the field coefficients c and rf, and another, exc{n{r))n{r) , is 
only known in terms of tabulated values. Expanding the product of two wavelets in terms of 
wavelets on finer levels would make possible the treatment of the cubic coupling to some level 
of approximation. (See, for example, 0.) However, this route becomes increasingly difficult 
for higher order interactions and is hopeless for non-algebraic or tabulated interactions, such 
as exc{n{r)). For higher order interactions it is natural, and for non-algebraic and tabulated 
interactions necessary, to evaluate the interactions at some set of points in space and then 
recover expansion coefficients for the result. One then relies upon the basis set to provide 
interpolation for the behavior at non-sample points. 

The benefits of both truncated wavelet bases and interpolation on dyadically refined grids 
are given by the use of interpolating scaling functions [1^, 0,0, [|14[) [@ (or interpolets 
P^ , [^), which are functions with the following properties (from pp. 6-7). 

Let 0(x) be an interpolet, then 

(INTl) cardinality: (p{k) = So,k for all e 

(INT2) two-scale relation: 0(x/2) = J2y<^z" Cy(j){x — y) 
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(INT3) polynomial span: For some integer m > 0, any polynomial p{x) of degree m 
can be represented as a formal sum X^yeZ" (^{y)4>{^ ~ v)- 

Cardinality allows the fast conversion between uniform samples and interpolating scaling 
functions and has subtle yet profound consequences for the resulting multiresolution basis. 
In particular, as is evident from our algorithms below, the expansion coefficient for a ba- 
sis function on a particular scale is independent of the samples of the function for points 
associated with finer scales. Consequently, the expansion coefficients which we obtain for 
functions maintained in our basis are identical to what would be obtained were function sam- 
ples available on a complete grid of arbitrarily fine resolution. This eliminates all error in 
the evaluation of non-linear, non-algebraic and tabulated interactions beyond the expansion 
of the result in terms of a finite set of basis functions. 

The Cy in the two-scale relation are referred to as scaling coefficients, and cardinality 
actually implies that Cy = (j){y/2). The two-scale relation allows the resolution to vary 
locally in a mathematically consistent manner. 

The polynomial span condition captures, in some sense, the accuracy of the approxi- 
mation. By cardinality, we actually have a{y) = p{y). We shall call m the polynomial 
order. 

Interpolets thus can be thought of as a bridge between computations with samples on 
dyadically refined grids and computations in a multiresolution analysis. The former point 
of view is useful for performing local nonlinear operations, while the latter is useful for the 
application of local linear operators. 

This manuscript explores 0{N) algorithms that calculate transforms and linear operators 
for grids of variable resolution but return, for the coefficients considered, exactly the same 
results as would be obtained using a full, uniform grid at the highest resolution without the 
need to introduce artificial temporary augmentation points to the grid during processing. 
We thus show that with relatively mild conditions on the variability of the resolution pro- 
vided by the grid, interpolet bases provided the ultimate economy in the introduction of grid 
points: only as many samples in space need be considered as functions used in the basis. The 
four transforms (forward, inverse, and the dual to each) mapping between coefficients and 
functions samples which we discuss here are particular to interpolet bases. For the applica- 
tion of operators in such bases, we show that the familiar non-standard multiply of Beylkin, 
Coifman and Rokhlin|^ shares with the transforms the property of correctness without the 
need to introduce additional grid points. Furthermore, we weaken the condition on grid vari- 
ability by using a modification of the non-standard multiply. We generalize the non-standard 
multiply so that communication may proceed between nearby but non- adjacent levels and 
thereby obtain less stringent conditions on the variability of the grid. All of theoretical re- 
sults in this manuscript are presented in a general d-dimensional space. Illustrative examples 
for the purpose of discussion will be given in d = 1 and d = 2 dimensions. The examples 
of applications in the final section will be in (i = 3 dimensions. Our focus is entirely on 
interpolet bases, and so it remains an open question whether these results hold true or can 
be adapted to more general wavelet systems. 

Our organization is as follows. In Sections 2 and 3, we explain how to manipulate and 
construct interpolet expansions and some aspects of how well they perform. These sections 
will present nothing new to the experienced wavelet user, but will explain our notational 
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conventions and recapitulates common theorems (|]T2[, [2^|, [13|, et al.) for wavelet 
novices. In Section 4, we describe how nonuniform bases can be conceptualized in the 
framework of interpolet expansions and then use our results to develop algorithms for the 
transforms. Section 5 details the algorithm for and other operators. Section 6 gives 
some practical details for the reader interested in implementing these algorithms. Finally, 
Section 7 compares, in three dimensions, timings of these implementations with the timings 
of naive algorithms. This final section also explores the convergence of a preconditioned con- 
jugate gradients algorithm in the solution of Poisson's equation for the full three dimensional 
electrostatic potential arising from the nuclei in the nitrogen molecule. 



2 Introduction to Interpolets 

There is a unique set of interpolets on R having symmetry and minimal support for a 
given polynomial order m = 2/ — 1 (the Deslauriers-Duhuc functions [|14|)- These are the 
functions with which this article is primarily concerned (our results carry over to more general 
interpolets, and no use will actually be made of minimal support or symmetry). 

To determine the Cj^'s, one sets C2j = Smo and Cy = c^y. One may solve the Vandermonde 
system, 



/I 
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• (2/ -1)2 
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to obtain the remaining Cy's. These coefficients satisfy the conditions for polynomial order 
21 - 1. 

The scaling coefficients for m = 1 are 



c_i = ci = 0.5, Co = 1, 



and for m = 3 (the example used for Figure 1.) they are 

_ _ 1 _ _ 9 _i 

C-3 — C3 — — — , C_i — Ci — — , Co — 1. 

Id Id 

One may then take tensor products of 0's and c^'s to form interpolets in higher dimen- 
sions. 



2.1 Interpolet Mult iresolut ion Analysis 

We are concerned with recursive representations of functions from samples at integer points 
on both uniform and refined grids. There are many definitions which make the exposition 
more clear. 

Definition 2.1 For k>0, let Ck = 2^2"^, and let Dk = Ck-i-Ck- For k<0, let Ck = 
and Dk = 0. 
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We consider to be the set of coarse lattice points on the lattice 2^~^Z'^ and D^, the 
detail lattice points, to be those points on 2'^~^Z" which are not coarse. Note: D]^ U — 
Cfe_i, and = CkU Df-U Dk_i U ■ ■ ■ U Di is a partition of Z". 

Definition 2.2 We let 9k{y) — min{k, m) where m is the largest integer such that 2^ divides 
all of the components of y. We call 9k{y) the level of the point y. 

Given the partition Z"" = CkU DkU Dk-i U • 

Ok{y) = I \[ 

Definition 2.3 Let S C Z". Let (f){x) be an interpolet. Letlk{(l), S) be the space of functions 
given by formal sums of the form J2v&s a(y)0( )• 

Where (p and S are understood, we may simply write X^. 

Definition 2.4 Let S C Z". Let J-k{S) be the vector space of R- or C- valued functions on 
which are zero at any point not in S (i.e. with support contained in S). 

Where S is understood, we may simply write JF^. Note: J^k{S) = Tk{S fl Ck) © ^k{S n 
Dk) © J'kiS n Dk-i) © • • • © Tk{S n L>i). 

The meaning of the k subscript will be established by the next definition, which will link 
vectors in JF^, with functions in X^. It is for this reason that while the jF's are technically 
identical, they are semantically different. In practice, the .F's are the actual data structures 
being stored on the computer. 

Definition 2.5 Let be an interpolet. Let if. : S ^ Ik{<l>, S) be defined by 

This definition extends linearly to the mapping <.| : J^k{S) — >■ Xk{(t>-iS) defined by: 

yes 

i.e. 

{itv){x)^ Y: v{ymx-y)/2^)-r E viy^x -y) 12^-') ■ ■ ■ E v{y)<i>{x-y). 
yeSnCk yeSnDk yeSnDi 



■ ■ U Di, we have 



y^Ck 
-1, ye A. 
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The set S can be thought of as the set of points in a refined grid. The identifications 
allow one to think of S* as a set of functions, {{L.f.y)\y G S}, which form a basis of X^.((/), S). We 
will sometimes refer to 5 as a refined grid and sometimes as a basis with this identification 
understood. 

One should think of the J^k as spaces of coefficients for function expansions in the corre- 
sponding Ik spaces, in the basis S. The if. simply associate a set of coefficients in J^^ with 
a function in X^. When is understood, we may write just l^. 

We are now in a position to state the basic theorems of interpolet expansions on uniform 
grids. 

Theorem 2.6 Let (f){x) be an interpolet on R^. Then each mapping : J-k{S) — > 1k{(t>-i S) 
(A; = 1, 2, . . .) is an isomorphism. 

Proof: Since the map ik is surjective, it is only necessary to show that ik is injective. 
By the definition, ikV = if and only if there exist v e J-k such that 

0= E v{y)4>{{x-y)/2^)+ ^ v{y)<P{{x - y) /2^-') + ■ ■ ■ + ^ v{y)<P{x~y). 
yeSnCk yeSnOk yeSnDi 

Let z e SnCk- By (INTl), we have (f){{z-y)/2'') = 5(2_y)/2fc,o = 3z,y, for y e SnCk, and 
also 4){{z-y)/2^) = 0, for y e SCiDi. So, {ikv){z) = v{z),z e SnCk, therefore v{z) = 0, z e 
S n Ck- This being so, one then has {Lkv){z) = v{z), z & S H Dk, so v{z) = z E S C] D^. 
Once again, (tfcf)(z) = v{z),z E S H -Dfe-i, thus we must have v{z) — 0, z E S H Dk-i- 
Continuing in this manner, we deduce that v{y) = Q,y E S, thus v = 0. 

□ 

Corollary 2.7 

ikis) = Xk{s n Ck) © Tk{s nDk)®-- -^Ms n d{) 

Since the sum is direct, the expansion is unique. 

This corollary is a consequence of observations in the above proof. 

Theorem 2.8 Let 4>{x) be an interpolet on EJ^ . Then 

yk,ik{<f>,z^)^ik-i{<f>,z^). 

Consequently, 

^h,k2,ikM z^)^ikM z^). 

Proof: To prove that Tk C Tk-i we note that Tk C Tk-i ^^k{Ck). Thus we just need 
to show Ik{Ck) C Ik-i- 

Translating by z E Ck and inserting powers of 2 where appropriate, one can rewrite 
(INT2) for (f) as 

<p{{x - z)/2'') = cpiix - z)/2'^-') + E Cy/^-^cpiix - z - y)/2'-'). 

yeDk 
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The terms in the right hand side are elements T^^i. Thus C T^^i. 

To prove C Xjt note that any element of Tj^_i can be expressed as: 

f{x)= E a{ymx-y)/2^-')+ ^ a{ymx - y) /2^-^) + ■ ■ ■ <y)^{^ ' v) 
yeCk-i yeDk-1 yeDi 

All the terms in this expansion but the first are elements of Ik- Since Ck-i — C^yj we 
may split the first sum up as, 

E a{ymx-y)/2^-') = ^ a(|/)0((a; - |/)/2^-^) + ^ aiy^x - y)/2'^-') 
j/eCfc_i yeCk yeDk 

The second term is also an element of Ik- 
Rewriting (INT2) one has {y e Cfc): 

- y)/2'-') = (piix - y)/2') - ^ c,/2.-i0((x - z - y)/2'-'). 

zeDk 

y e Ck, z e Dk, so y + z & Dk, thus the right hand side is made up of elements of Ik- Thus, 

Ik-i C Ik- 
□ 



2.2 Interpolet Transforms 

Corollary 2.9 Interpolet Decomposition 

The set of isomorphisms ik induces a set of isomorphisms 

Jki,k2 • 'F ki ^ ki,i 

JkiM = ^k2 ° ■ 

We refer to these isomorphisms as interpolet transforms. It is our convention to let Jk = Jo,k 
and J_fc = Jfc 0- The reader will note that the Jj are linear transformations on the coefficient 
spaces, and are thus the primary object of computation. 

We now turn to a study of the J's. It is clear from the definition that for ki < 
JkiM = Jk2-i,k2 o Jk2-2,k2-i o • • • o JkiM+i, and similarly for ki > k2- Thus, we need only 
study the Jfe,jfc+i and Jfe+i,fe mappings. 

Theorem 2.10 Computation Theorem 
Letv e TkiZ"")- 



{Jk,k+iv){y) = 

where v'{y) = v{y) - 'ZzeCu+^ C{y-z)/2''V{z)- 



v{y), y ^ Dk+i 
v'{y), y e Dk+i 



I k+i,k lyyi I ^,^y^^^ y g ^^^^ 



where v'(y) ^v(y) + EzeCk+i C{y+z)/2kv(z)- 



Proof: For v e we have 



expanding the first term, 

tkV= E v{y)cl>{{x-y)/2')+ J2 v{y)cl>{{x - y)/ 2') + J2 vivMi^ - y)/^''') + ■ 

using (INT2), 4>{{x - z)/2^) = 4>{{x - z)/2'=+i) - E,eD,+i c^y-z)/i^4>{{x - y)/2^), 

ikv= E v{ymx-y)/2^+')+ v'iy^x - y) /2^) + viy^x - y) /2'-') + 



{t'k+if'kv){y) 



v{y), y ^ Dk+i 
v'{y), ye Dk+i 

where v\y) = v{y) - EzeCk+i C{y-zy2^v{z). 

The proof for Jk+i,k is similar. 

□ 

Similar to what one might get with wavelets, we see that we can compute the coefficients 
of intcrpolet expansions on uniform lattices by a pyramid algorithm. Computationally, 
this procedure can be carried out by first computing the Di coefficients with Jo,i, then by 
computing the D2 coefficients from the Ci data with Ji 2, and so on. In this sense, it is no 
different from standard multiresolution decompositions. 

A feature of the interpolct decomposition is that the transformations all have a particular 
lower triangular form. That is, if we write v G J-'k as a vector with its Ck+i components first, 
its D^+i components second, and the rest of its components third, then the transformation 
takes the form, 

/ 7 0\ /vc,^^\ 
M / vd, 
/ vd, 
/0 
\ // \ . 

The inverse, Jk+i,k, is obtained by replacing M with —M. 



Jk,k+lV 



k + l 



3 Accuracy of Interpolet Approximation 

Given a function f{x) on /?", one can form an interpolet approximation to / by the formula: 

fix) ~ E f(jj)(Pix -y) = ^0/, 

yeCo 

where / on the right hand side is thought of as a function restricted to Co = (a more cum- 
bersome but more precise notation is io{/(z/)}|2/6Z")- This approximation has the property 
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Starting from the expansion 60/ G Xo(0)-Z^") one can construct equivalent expansions, 
i^kiJkf) £ 2rfc(0, Z"). The coefficients J^/ are referred to as the interpolet transform of the 
function /. 

If f{x) is sufficiently smooth, then we can expect that the coefficients {Jkf){y),y G -D;, / < 
k, will be small. This statement is captured rigorously by the following lemma and theorem. 

Lemma 3.1 Let 4> be an interpolet with polynomial order of m then 

p{x) e lN{(f>,CN) 

for any integer, N, and any polynomial, p, of degree m. 

Proof: p{2^x) is a polynomial of degree m. By (INT3), p{x) can thus be represented by a 
formal sum in To(0, Co), namely 

p(2^a:)= ^p(2^|/)0(a;-y). 

By changing variables, we may rewrite this as 

p(x) = p(2^y)0(x/2^ - y) 

= E P(2^2/)'/'((^-2^2/)/2^) 
= j:p{ymx-y)/2'') 

V&Cn 

yeCN 



□ 



Theorem 3.2 (Residual Theorem) 

Let f{x) be a polynomial function of degree m. Let (f) be an interpolet with a polynomial 
order of m. Then 

' f{y), y G Ck 
0, yeDi 



iJkf){y) 



Proof: 

By (INT3), f{x) G Xo(0, Co). By the lemma, we also have f{x) G Tk{(f),Ck). Recalling 
that Jkf gives the unique expansion coefficients of f{x) in the decomposition of Xo(0, Co) 
given by Tk{4>, Ck) Q)Tk{4>, D^) ® ■ ■ ■ Q)Tk{4>, -Di), we see that the Xfe(0, Di) coefficients must 
vanish, while the Tki<P, Ck) coefficients are given by the lemma, namely f{y) for y G Ck- 

□ 

The coefficients {Jkf){.y)iy G Di^l < k, are called residuals. The Residual Theorem 
suggests that the magnitude of the residual coefficients (Jkf) at a point y E Di are expected 
to be proportional to the (m + l)th derivative of f{x) at y. See Figure 2. 
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4 Truncated Bases 



Typically, one truncates an expansion by eliminating elements of a basis, setting their coef- 
ficients to 0. One then said to be working in a truncated basis when one works within the 
subspace formed by the remain basis elements. In the notation of this paper, this corresponds 
to taking expansions in Tkid). S) with coefficients in Tk{S). 

One may also view a truncated basis as the set of expansions one gets when the coefficients 
of some of the basis elements have been set to "don't care" values. Mathematically, this is 
accomplished by quotienting out the "don't care" elements. 

Definition 4.1 Let 

Xf =X,(0,Z-)/X,(0,Z--5) 

and 

When the identification of F^^iS') with can be made will be dealt with later in this paper. 
For now, one may view these definitions as a trick to make the proofs less complicated and 
for understanding exactly why and in what sense the algorithms are correct. Once again, we 
think of JFf as a grid on which the elements outside of S have "don't care" values, and Tk{S) 
as a grid on which the elements outside of S vanish. The continue to be isomorphisms 
(since ik{J'k{Z'' - S)) = 2:^(0, - S)). 

However, it is not necessarily true that X^^ — X^^. When this condition fails, then it is 



-1 



2 fe2 



no longer possible to define Jk-^ 

To be sure, one could still define some sort of Jki,k2 by setting the elements in Z" — 5" 
to zero, then applying the full grid version of Jki^i ^^'^ then considering only the elements 
in S of the answer. This definition by itself has some drawbacks. Mathematically speaking, 
this is the same as JkiM = P ° l-^I ° ^ki ° r, where r : J^f^ — > J^ki, is some lift to the full 
grid, and p : X^^ Xf^ is the standard projection operator onto the quotient. Generally, if 
one follows this approach one will no longer have Jki.k-z = Jk^-iM ° Jk2-2,k2-^ o • • • o Jki,ki+i 
because Xf are not all equal. In terms of diagrams, where we once had 

f T — T T 
fci -^k\ -^k2 *' fc2 5 

we now have 

-r-S ''H q-S J_ q-S -r-5 
•^fei ^ -^fei — -^ki ~^ ki- 

What one needs is a condition on S such that Xk^iZ'^ — S) — Xk^iZ"^ ~ S). If this were 
true, then the definition of the operator as Jki,k2 = ^'kl ° '-fci would actually be independent 
of the values of the elements of Z" — S. In that case the quotient spaces are identical. 

Definition 4.2 We say that the set S is a good basis when it satisfies the condition \/ki, k^iXk^ {Z^- 
S) = Xfe,(Z" - S), and thus Xf^ = Xf^. 
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To get a handle on this definition, one sees that this is achieved when Tk{Z^ — S) = 
Tfc+i(Z" — S). For this to be so, whenever y ^ — S, every z such that 0(x/2^^^^^ — z) is 
in the two-scale expansion for 0(a;/2^'^^^) — y), must also be a member of Z" — S. 

This can be captured in the following table (in which we let 9k{y) = 9k{z) + 1). 



in expansion? 


zeS 


zeZ'^-S 


yes 


ok 


ok 


yeZ'^-S 


not ok 


ok 



The good basis condition for fast synthesis and reconstruction has also been discovered by 
Cohen and Danchin (see S-trees in a coming work||10||) which appeared after the original 
submission of our manuscript. 

For some of the algorithms presented, we may employ additional conditions based on the 
supports of the functions themselves (not just the support of their expansions). 

Definition 4.3 We say that functions f and g touch whenever supp{f} fl supp{g} has 
nonzero measure. 

For p > we say that S has the p-level touching property when it satisfies the condition, 
that for y E Z'"' — S, z E 2'", 9k{z) < 6k{y) — P, o-nd L^y touches l^z implies z & Z^ — S. 

A less formal way of phrasing this definition for the case of 1-level touching is that if a 
level / point,?/, is a member of — S then any point, z, at level / — 1 or lower for which tkU 
touches ikZ must also be in — S. For 2-level touching, one only considers any points at 
level / — 2 or lower, and so on for p-level touching. 

The allowed touching possibilities can be summarized the following table (in which we 
let ek{y) > Okiz) +p). 



touch? 


zeS 


z e - 5 


yes 


ok 


ok 


yeZ'^-S 


not ok 


ok 



One example of a 1-level touching good basis for one dimensional 3rd order interpolets 
(the ones being used as an example) is the set of interpolets centered at the points 

/n-l 

S = CnU i\J{-7 ■2\-5-2\...,5-2\7 ■ 2'} 

\/=o 

The support of the interpolet on level at y is [—3 + y,3 + y], the union of all the supports 
of the interpolets in 5* on level is [—10, 10]. The support of an interpolet on level 1 at |/ is 
[—6 + y,Q + y], thus the only interpolets on level 1 which touch the interpolets on level are 
precisely those ones at points — 14, . . . , 14, which are precisely the ones included in S. No 
interpolet not included on level 1 touches an interpolet included on level so the definition 
is satisfied. The argument proceeds similarly on higher levels. In three dimensions, this 
example corresponds to nested concentric cubes of size 15^ ■ ■ • 2' at each level / < n. 
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An example of a 2-level touching good basis for 3rd order is the set of interpolets centered 
at the points 

S^Cnu(\j{-5-2\-3-2\...,3-2\5-2'}j . 

Note (for n > 2) that this set is not 1-lcvcl touching since the level 1 interpolct centered at 
7/ = 12 is not included, while an interpolet it touches, namely the level interpolet centered 
at y = 5, is included. 
Also, the set of points 



^ = C„U (|"U{-3-2',-l-2',l-2',3-2'}^ . 



forms a good basis, but is not 1-level touching or 2-level touching (it is 3-level touching, 
though) . 

The above examples are meant to suggest that the good basis and the p-level touching 
definitions can be thought of, informally, as conditions telling one how quickly one can change 
from one resolution to another. Essentially, any nested set of refined regions can satisfy these 
conditions so long as the margins around the set of points at a given resolution are wide 
enough. 

^From a computational point of view, what these conditions do is ensure that data-paths 
which carry coefficient information between different resolutions are not broken by zeroed 
coefficients at intermediate levels. 

It is clear from the preceding discussion, in a good basis, one has JkiM — JkQ-iM ° 
Jk2~2,k2-i ° • • • ° -^fci,iki+i- We may now generalize the computation theorem to a truncated 
basis. 

Theorem 4.4 Good Basis Computation Theorem 

Let S be a good basis, v e !Fk{S), y E S 

and V is any member of the equivalence class of v 



{Jk,k+iv){y) = 
where v'{y) = v{y) - Ezesnc.+i c^y~z)/2M^) 

iJk+i,kv)iy) = 
where v'{y) = v{y) + EzeSnCk+i C{y+z)/2''V{z) 



v{y), y e S - Dk+i 
v'{y), y e Dk+i 



v{y), y e S - Dk+i 
v'{y), y e Dk+i 



Proof: 

In a good basis, the computations of all Jki,k2 ^i"^ independent of the representative. 
Thus, this algorithm computed on v gives a member of the same class as would be computed 
on V. 

□ 
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Thus, the pyramid algorithm of the uniform basis, Z", has a counterpart in a good basis 
S, allowing the computation of the expansion coefficients in J-'k{S) from the values of the 
function in J-'o{S) (and also has the lower triangular structure). 

In a good basis, one thus has the ability to perfectly reconstruct the multiscale coefficients 
of a function for the basis functions associated with the points of the refined grid 5" by simply 
applying the pyramid algorithm on zero-lifts at each stage of the algorithm. The above 
theorem establishes this as true, even though we do not necessarily expect the data zeroed 
during the lift to be small. (The function may have significant sample values throughout 
the domain of the representation). Also, with exact recovery of sample values, it is easy to 
perform local nonlinear point-wise operations of the form f{x) G{f{x)) (e.g. e^^^^), or 
point-wise multiplication (i.e. f{x),g{x) f{x)g{x)), despite the truncation. 

The reader may note that this result is " analysis- free" in that we have sparsified the com- 
putation, not by proving that the coefficients vanish outside of the truncation for some class 
of functions, but by showing the the coefficients we wish to compute have no dependency on 
the coefficients we omitted. Computationally, this means the data-structure in the computer 
requires no intermediate augmentations (contrast with |TB[). 

The advantages conferred by the additional the p-level properties are seen in the context 
of operator evaluation, and will be the subject of the next two sections. 



5 Multilevel Algorithms for and Other Operators 

Given v,w E J^k{Z"-), we may compute the stiffness matrix of a model operator, V^, 

(^LkV\V'^\Lkw'^ = J {ikV){x)\/'^{ikW){x)(r-X 

by changing the expansions, 

(^LoJ-kv\V'^\ioJ-kw'^ = {J~kv){y){J-kw){z) (f){x ~ y)V'^(f){x - z)d''x. 

This reduces the computation to the computation of the matrix elements {4>{x — {y — z)) | V^|</)(x)) 
which can be done by solving the associated eigen-problem obtained by applying (INT2). 
In particular, let L° = {(f){x — ?/)|V^|0(x)), then 

Ll= Y c,,c,,(0(2x-2y-^i + ^2)|V2|0(2x)) 

rO _ 92-n rO 

zi,z2eZ" 

which we solve by standard methods (found in [^, for example). In subsequent sections 
of this article, we will also define L+ = (0(x - y)\V'^\4>{x/2)), Ll^^ = (0(x - ?/)| V^|0(a;/4)), 
Ly = (0((x — ?/)/2)| V^|0(x)), and = (0((x — ?/)/4)| V^|0(x)), which can be computed 
from L° by employing (INT2). Although it is true by Hermiticity that L~ = Ly and 
L~~ = Ly~^, we will make no use of this fact. 
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One can write a matrix expression {tkv\V^\Lkw) — v^Jt^LJ-kW, where the L is the 
Toephtz matrix with coefficients Ly. In practice, one typically formulates the preceding as 
a computation oi u = J^_j^DJ_kW for some w G Then {Lkv\W'^\Lkw)= v^u. u may be 

thought of as an element of J-'k{Z"-)*, the dual space of J-'k{Z"'). The task is then to compute 
the coefficients u{y),y e Z'" — Ck®Dk®- ■ - ^Di. Any algorithm for computing x*Ay can be 
adapted to an algorithm for Ay, and for purposes of making proofs, it is somewhat easier to 
keep thinking of the computation as {Lkv\'V'^\ijsw), which is the point of view we shall take. 

However, computing {Lkv\'V'^\L.kw)= v^J'ij.LJ^kW, by just applying the transforms, and 
the Toeplitz matrix is problematic, since this process makes it necessary to either represent v 
and w on a uniform grid, or to compute a matrix element between each pair of functions in the 
truncated expansion which touch. In the first case, one ends up with an 0{N) computation 
for (typically) a very large N. In the second case, one chooses between extremes which are 
0{N) and quite complicated or simple and 0{N^). 

The following sections outline the design of multilevel algorithms for for both 1-level 
and 2-level touching bases. Both algorithms are derived according to the following format: 

break up the expansion of {Lkv\V'^\Lkw) into a decomposition over elements at the same 
level and adjacent levels. 

rewrite the expansions in terms of the matrix elements between elements of those levels 
and the transforms of higher/lower level elements. 

implement the algorithm by computing those terms separately. 

estabhsh correctness in a p-level truncated basis. 

Although only the 1-level and 2-lcvel algorithms have been explored in any detail, this 
same process will generally work to produce 0{N) p-level algorithms for any p. 



5.1 in 1-level decomposition 

The 1-level decomposition of {Lkv\V'^\Lkw) is 



^ikv\V likwj = (^ikv\V likwj + (^ikv\V likwj + (^ikv\V likwj , 
where 

Lkv\V^\ikwy = {iky\^\kz) v{y)w{z) 

ikv\V'^\ikwy = Yl {i^ky\^'^Vkz) v{y)w{z) 
Ok{y)<Ok[z) 

ikv\V'^\ikw) = Yl {^ky\^'^Vkz)v{y)w{z). 

That is, we express the product as contributions from levels to the same level, higher levels, 
and lower levels. We will now investigate each of these terms individually. In the language 
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of matrices, these respectively correspond to diagonal, above diagonal, and below diagonal 
blocks of the matrix. 



l'kV\V'^\l'kw) = f* 



({\f 


(1) 


(1)" 


(1) 


(1)" 


(1) 


\{\r 


(1) 



/i\0 



w 



5.1.1 Diagonal blocks: (l)'^ 

For some fixed 1,0 < I < k a. term of the form 

E - y)/2O|V^|0((x - z)/2^)) v{y)w{z) 

dk{z)=9k{y)=l 

contributes to (|)°. 

However, - y)/2^)\V^\(P{{x - z)/2^)) = 2'("-2) {4>{x - y)\V'^\(t){x - z)) = 2'("-2)L0_^. 

Thus we have 



5.1.2 Super-diagonal blocks: (|)''" 

For some fixed 1,0 < I < k a term of the form 

Y: {Hi^ - y)/2^)\V'\cl>{{x - z)/2^')) v{y)w{z) 
9k{z)=l'>l=9k{y) 

contributes to (l)"^. 

Applying the inverse transform ( Jj. ;+i) to the w coefficients in the sum allows us to write 
this term as 

Y: {<t>{{x - y)/20|VX(x - z)/2^+')) v{y){J,,^,w){z). 

6i+i{z)=l+l>l=6k{y) 



Thus we have 



z . 



5.1.3 Sub-diagonal blocks: (|) 

For some fixed 1,0 < I < k a term of the form 



E - y)/2^')\V'\<P{{x - z)/2^)) v{y)w{z) 

ek(z)=l<V=9Uy) 



contributes to (|)" 
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This time applying the inverse transform (J^ ^+i) to the v coefficients in the sum allows 
us to write this term as 

E {<l>{{x - y)/2^^')\V'mx - z)/2')) {J,,^,v){y)w{z). 

9k(z)=l<l+l=9i+i{y) 

Thus we have 

{i,v\^%wy ^Y.'^'^--"^Y.L:-y{Jk,i^iv){y)w{z). 
5.1.4 Implementations 

The above observations demonstrate the correctness of the following algorithm: 

input: v,w & !Fk{Z'^) 

output: ans = (6fcv|V^|tfcw) e R 

Let wtmp = w, let vtmp — v, let ans — 

for / = to A; 

ans ^ ans + 2'(-2) E^^(^)=^^(,)=^ Ll_yV{y)w{z) 

end for 

for Z = A; — 1 down-to 

ans ^ ans + 2'("-2) Lt_yV{y)wtmp{z) 

wtmp <— Ji+i^i{wtmp) 

end for 

ior I — k — 1 down-to 

ans ^ ans + 2'("-2) Ee,^,(y)=i+i,e,{z)=i L-_yVtmp{y)w{z) 
vtmp <— JiJ^i^i{vtmp) 

end for 

Note, we have made use of the fact that Jk^i = Ji+i,i ■ ■ ■ Jk,k-i so that at the beginning of 
each iteration in the second (last) loop wtmp = Jk^i+iw {vtmp = Jk,i+iv). 
We observe that this algorithm is 0{N) in time and space. 

We may adapt this to an algorithm to compute u e ^^(Z")* such that u{v) — V^lifety) . 
input: w eJ^k{Z'^) 

output: u such that u{v) = {tkv\V'^\Lkw) 

Let wtmp = w, let M = e J^kiZ"")*, let utmp = e ToiZ"")* 
ior I — to k 

u{y) ^ u{y) + 2'("-^) E^, Ll_^w{z) 

end for 

for I = k — 1 down-to 

u{y) ^ u{y) + 2'("-2) i:e,^y)=i,ievek+,^z)=i+i Lt-yWtmp{z) 
wtmp ^ Ji^i^iwtmp 

end for 

for Z = to A; - 1 

utmp Jf^^ jUtmp 

utmp{y) ^ utmp{y) + 2'(»-2) E^,^^(^)=;+i,^^(^)=i L-_yW{z) 
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end for 

■u M + utmp 

The third loop is the result of transposing the linear operator Ji+i^i in the last loop in 
the previous algorithm. We have also used the fact that J\i = ■ ■ ■ Jf^^ i to ensure that 

at the beginning of each iteration in the third loop, utmp G J-'i+i{Z"')* . 

It is easy to check that this is an 0{N) algorithm in both time and space. It is very 
important for atomic structures computation that this algorithm scales linearly with the 
number of atoms. Without such a scaling, one can only compute electronic configurations 
for small molecules. 

The reader may note a similarity between this algorithm and other matrix-vector mul- 
tiplies used to apply operators in a uniform wavelet basis. In fact, the 1-level algorithm 
presented above is identical to the nonstandard multiply found in Q and developed for or- 
thonormal wavelet bases. The nonstandard multiply was introduced by Beylkin, Coifman, 
and Rokhlin to sparsify integral operators whose kernels were smooth or vanishing off the 
diagonal, while keeping a uniform basis. 

However, in contrast to that program of sparsification, interpolets allow one to sparsify the 
basis, and, with out introducing additional grid points, still be able to apply the nonstandard 
multiply routines, with any local operator. With interpolets, we remove any elements from 
the expansion that we believe will be insignificant, still having a good approximation to 
our function at the points we retain. Beylkin et al. express the matrix elements of the 
operator itself in a nonstandard orthonormal basis and then remove those matrix elements 
which are determined to be very small to produce a sparse matrix. 

The use of interpolating scaling functions has achieved some degree of simplicity and con- 
venience in carrying out fast point-wise operations. Although there is no associated difficulty 
in electronic structure calculations [|I|, for other applications, the loss of orthogonality might 
be too great an expense. In those cases, one might consider employing compactly supported 
approximations to orthogonal interpolating functions found in |]^ . It appears that with some 
additional complexity one might be able to extend the present algorithms to other wavelet 
bases. The additional complexity of other schemes and the need for fast point-wise opera- 
tions in our applications are the chief reasons we do not consider doing this in the present 



work. Finally, there is a large body of work the reader may wish to consult (0, ||T^, [|TT 



21| , and for adaptive refinement techniques when, in contrast to the case of electronic 



structure calculations, the behavior of the needed refinement is not known a priori. 



5.1.5 Correctness in a 1-leveI touching good basis 

The above decomposition of the product and the associated algorithm is what we seek to 
extend to a good truncated basis. In practice, one takes the zero-lift representatives of v 
and w G JF^ and computes {Lkv\S/'^\Lkw) . By the computation theorem of good truncated 
bases, the value of {Jki,k2'^)iy),y ^ is independent of the representative (likewise for w), 
however, we must also address the issue that {Jki,k2'^){y) 7^ 0,?/ ^ S' (i.e. JkiM''^ JkiM'^)^ 
and thus y ^ S may have a contribution to the decomposition above, requiring us to augment 
S in order to get the right answer. 

Theorem 5.1 If one replaces with S everywhere in the Multilevel algorithm for {Lkv\'V'^\Lkw) , 
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then the algorithm computes 

Proof: As mentioned in the remarks, the multilevel algorithm requires Jk,i = Ji+i,i ■ ■ ■ Jk,k-i, 
which is true in a good bases. 

The (|)° computation proceeds identically for either Z"' or S, so the first loop will con- 
tribute correctly the term V^jifctt')''. 

To check the (|)^ contribution, one observes that the necessary term is 

{ikvlV'hwy = J:2^^^-'^ Y: L-;_^{.h,+iv){y)w{z). 

' yeZ'^,zeZ'^:9l+i{y)=l+l,9k{z)=l 

Suppose that 3y ^ S and z E S such that L^^^ ^ 0. This implies that supp{<./+i|/} fl 
supp{tfc2;} has nonzero measure. Since supp{ii+iy} C supp{tfc|/} we conclude that supp{<.fc|/}n 
supp{tfe2;} has nonzero measure. This cannot be so if 5" is 1-level touching, and since 
^(z) = 0, z ^ 5" we may restrict the sums over y and z in the contribution, 

2'^"-'^ E L-_,{Ju,^,v){y)w{z). 

y€S,z€S:9,+i{y)=l+l,9k{z)=l 

The proof for (|)~ is identical with v's and w's reversed. 

□ 

Immediately we have the following: 

Corollary 5.2 // one replaces with S everywhere in the Multilevel algorithm for u such 
that u{v) = {Lkv\V^\Lkw) , then the algorithm computes 

The computation of {Lkv\V'^\Lkw) serves as a template for another common computation 
one may wish to perform, namely = j{ikv){x){ikw){x)dr'x, i.e. the L'^{BJ^) inner 

product of tfeV and i}^w. 

5.2 Computing Other Operators 

To compute {ikv\ikw), one simply replaces L°, L"*", and L" with = {(j){x — y)\(j){x)), Gy = 
{(t){x-y)\(t){x/2)), and Gy = - y)/2)|0(a:;)), and then replaces the factors of 2'^"-^) 

with 2^". After that, the algorithms and theorems for carry over directly. 

The above procedure can be used for creating a multilevel algorithm for any operator, 
C, which is 

local :supp{C/} C supp{/} 

translation invariant :Of{x + a) — {Of){x + a) 
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homogeneous :Of{sx) — s°'{Of){sx) 
by forming the appropriate coefficients, 0^^'^~^, and inserting appropriate factors of 2'^"+'^). 

However, locahty is the only property which is really required for 0{N) multilevel algo- 
rithms so long as one can compute O/^^^/ ^ efficiently. 

As an example of a local, homogeneous, but not translationally invariant operator, we 
shall discuss the coefficients for the multilevel algorithm for the xi operator in two dimen- 
sions. 

We first consider the coefficients Xij^^^^^ given by 
= / <l>{{x,-m^)/2^^+'^)<|>{{x,-m,)/2^^+'^)x,^^^^ 

i.e. 

x^n,m,m' = / - mi)/2')0((x2 - m2)/2')xi(j){{xi - m[)/2')(t)({x2 - m'^)/2')dxidx2 
^\i,m,m' = / 0((^i - mi)/2')0((x2 - m2)/2')xi0((xi - m\) /2'+^)(l>{{x2 - m'^)/2'+')dx,dx2 
£i7,m,m' = / 'Pii^i - mi)/2'+i)</.((x2 - m2)/2'+i)xi</.((xi - mi)/2')</.((x2 - m'^)/2^)dxidx2. 

Separating the xi and X2 integrations, we see from this that we may write afi/^^ ^ = 

(mi— ni)/2' i,m2,n2 



IX 

and G is defined above. The problem of computing the xi coefficients has been reduced to 
computing the X coefficients and then doing a multiply with the already known G coeffi- 
cients. 

In addition, one has 

J (piix- m) /2^'-+^^)x(l){{x - n) /2^^+^^)dx ^ n J (l){{x + n - m) /2'(+^))(/.(x/2'(+^))cix+ 

(j){{x + n- m)/2'(+^))x0(a;/2'(+^))rfa;. 



and thus 



+,-} 

)/2l 



1 4>{{x - m)/2'(+^))x0((x - n)/2'(+i))dx = n2^GfC.nW + '^"'^t- 
where 

Sy — j (j){x — y)x(j){x)dx 

Sy — J (p{x — y)x(f){x/2)dx 
= J cj>{{x-y)/2)xcj>{x)dx. 

Thus we see that for the xi operator, 

. {0,+ ,-} _ r,i^{0,+ ,-} /'„9/r'{0,+ -} , 92ic{0,+ -}n 

giving an efficient means to compute the multilevel coefficient for this operator. 
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5.3 in 2-level decomposition 

The previous algorithm for 1-level touching good bases can be expanded to 2-level touching 
good bases. One may wish to do this because one finds that the 1-level touching property 
is too stringent and requires one to augment one's basis set far too much to be practical 
comput at ionally. 

Much of the reasoning for the 2-level case can be found in the details of the 1-level case, 
so the exposition here will be more compact. The resulting algorithm will be correct for 
2-level touching good bases. 

The 2-level decomposition of {ikv\V'^\ikw) is 

+ {LkV\V'^\ikW)^^ + {ikV\V'^\ikW) 

where 



ikv\V \ikW 



ikv\V \ikW 



0M+i=9k{z) 



++ 



9k{y)-l=9k{z) 
9k{y)+l<9k{z) 
9k{y)-l>9k{z) 



Which corresponds to the matrix decomposition: 
ikvlS/'^kkW 







v{y)w{z) 






v{y)w{z) 




i^kz) 


v{y)w{z) 




^kZ^I 


v{y)w{z) 




I'kZ 1 


v{y)w{z) 



1 


(1)^ 




(1)^" 


\ 


{\r 




(1)" 






(1)" 


{\r 








V(l)" 


(i>" 


(1)" 




/ 



w. 



The key idea is to evaluate the diagonal and first off diagonal blocks of the matrix and 
then to compute the other blocks above and below the tridiagonal through the transforms. 

5.3.1 (|)°, (|) + , (I)-, (|)++, and (1)" 

The definitions for the contributions in the decomposition proceed just as they did for the 
1-level case. 



ikv\V\kw) =^2'('^-2) Ll_^v{y)w{z). 
' 9k{z)=9k{y)=i 
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' ek{z)-i=ek{y)=i 

I 0k{z)=0k{y)-l=l 

Lkv\V'hwy^ = ^2^("-2) Y: Lt^yV{y){Jk,i+2w){z). 
I 0k{z)-l>l=9k(y) 

ikv\V'\ikwy~ = ^2^("-2) Y L--y{Jk,i+2v){y)w{z). 
I eaz)=i<9k(y)--^ 

5.3.2 Implementations 

input: v,w & 

output: ans = {Lkv\V^\Lkw) G R 

Let wtmp = w, let vtmp — v, let ans — 
for / = to /c 

ans ^ ans + 2'(-2) E^^j,)=^^(,)=; L^,_yV{y)w{z) 

end for 

for Z = to A; - 1 

ans ^ ans + 2^(^-2) E^,(^)=;,^^(,)=;+i Lt_yV{y)w{z) 

end for 

for / = to A; - 1 

ans ^ ans + 2'("-2) L-_yV{y)w{z) 

end for 

for Z = A; — 2 down-to 

ans 4- ans + 2'("-2) ^^^^^^^^^^^^^^^^^^^^ L+4^(y)wimp(z) 
wtmp <— Ji+i^iwtmp 

end for 

for Z = A; — 2 down-to 

ans ^ ans + 2'("-2) E^^^^(^)=;+2,^,(.)=/ L-IyVtmp{y)w{z) 
vtmp <— Ji+i^ivtmp 

end for 

We adapt this algorithm to compute u e Tk{Z'^)* such that li(i') = (t^f | V^|tfcw) 
input: e 

output: u such that = {ikv\V'^\Lkw) 

Let tytmp = w, let M = e :^fe(Z")*, let wimp = e :Fo(Z")* 
for / = to A; 

u{y) ^ u{y) + 2'(-2) E^,(,)=^,(.)=, 

end for 

for Z = to A; - 1 

u{y) ^ u{y) + 2^(-2) E^,(,)=,^,(,)=,+i i^t.^(^) 
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end for 

for Z = to A; - 1 

u{y) ^ u{y) + 2'("-2) Ee,^y)=i^,Ai^)=i L-_^Hz) 

end for 

for Z = A; — 2 down-to 

u{y) ^ u{y) + 2'(-^) T.eay)=ifi,,,iz)=i^2 Lt^y^tmp{z) 
wtmp <— Ji^i^iwtmp 

end for 

for / = to A; - 2 

utmp ^ Jf^i lUtmp 

utmp{y) ^ utmp{y) + 2^^^-^^ ^e^+2{y)=i+2,e,iz)=i L',IyW{z) 

end for 

u ^ u-\- utmp 

6 Efficient Implementation 

We have produced a very successful 3D implementation of all of the above algorithms for the 
interpolet used as this paper's example (m = 3). Implementation details are given in this 
section. The ideas used to make this implementation efficient for all of the above algorithms. 

The purpose of this section is to give additional information to readers who wish to 
implement these algorithms themselves. 

6.1 Data Structures 

The interpolet data and function samples are kept in a sequence of blocks at various levels. 
Each block at level k contains the points of a rectangular subset of Ck-i- Since D^. = 
Ck-i — Ck, we use the collection of blocks at level k < p {p being the top level) to represent a 
rectangular subset Ok, ignoring the Ck points of each of these blocks. In our implementation, 
these extra Ck points hold the value in between operations and take on useful intermediate 
values during operations. Since we are working in 3 dimensions, this multiplies the storage 
required by a factor of about |, which we found an acceptable cost for its advantages. 

The coefficients for the transforms and operators are kept in various 3D arrays. Al- 
though it is possible to build the coefficients upon demand from a set of ID arrays of 
coefficients, we have found that the arithmetic cost of doing this is much greater than the 
cost of storing them (about 10 flops are required for each coefficient, while the 3D ar- 
rays are still small enough to be stored in cache). We have (in Fortran notation) the filters 
cs(0:3,0:3,0:3), SAMELEVEL(0:5,0:5,0:5), ONELEVEL(0:8,0:8,0:8), and (for 2-level algo- 
rithms) TWOLEVEL(0:14, 0:14, 0:14) (note: we have made use of the fact that our operators 
are symmetric to cut the size of these arrays by a factor of | and use ONELEVEL and 
TWOLEVEL for both upward and downward inter-level communication). 
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6.2 Implementation Overview 

The blocks described in the previous section are used as the fundamental objects for manip- 
ulation. The computation proceeds by employing block-to-block subroutines for the various 
operations, having every block at a level send data to every block at the same level, one level 
up or down, or (for 2-level algorithms) two levels up or down. 

The number of blocks at each level is not very large, and if a subroutine determines 
that the intersection of two blocks is empty (which it does by examining the bounding 
rectangles), then it returns immediately. Thus, while this algorithm is to be 0{B'^) where B 
is the number of blocks, it remains 0{N) where N is the number of actual points, because 
B is much smaller than N. 

6.3 block-to-block Subroutines 

The block-to-block subroutines arc all designed to take two blocks (source and destination) 
and a set of filter coefficients and place the result of convolving the filter with the source 
block in the overlapping points of the destination block. There is a block-to-block subroutine 
for the interpolet transform, its transpose, its inverse, and its inverse transpose, as well as 
operator application routines for the same-level operator, up-one-level operator, down-one- 
level operator, up-two-level operator, and the down-two-levcl operator. 

All of these routines precompute the bounding box for the points in the destination 
block which are in the range of influence of the source block and, for each point in this 
sub-block, the bounding box for the points in the source block in the domain of influence 
of the destination point. The result of this precomputation is that the only data values of 
the source (destination) which are accessed are the ones which are read (modifled). This 
decreases the number of data accesses in our test problems by a factor of 7. 

Additionally, blocking the computation generally increases the locality of access for the 
data. More data requests hit the cache than would occur in a more arbitrarily arranged 
construction. 

7 Results and Conclusions 

With interpolets, it is possible to carry out 0{N) computations in truncated bases, where 
N is the number of elements retained in the truncation, without having to augment the 
grid of points associated with the functions maintained in the basis. Along with allowing 
one to compute common linear physical operations, interpolet algorithms also allow one to 
transfer between function values and multiscale expansion coefficients on grids of variable 
resolution recovering the same results as one would obtain working with data on a full grid 
of arbitrary resolution but without introducing additional grid points into the calculation. 
This allows local nonlinear couplings to be computed quickly without the introduction of 
additional sample points and without the introduction of additional approximations which 
must be controlled. 

These algorithms have been implemented in FortranQO and have subsequently been 
adopted for use in electronic structure computations as described in the introduction. Prior 
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to this, we had been using very simple, naive 0{N'^) algorithms which implement each trans- 
form and operator as multiplication by the multiscale representation of the corresponding 
matrix. These multiplies check all points for being within interaction range and then con- 
struct the appropriate matrix element as needed. This is required in the naibe approach 
because the variety of inter-scale matrix elements is too wide to store in table of reasonable 
size. This algorithm ultimately scales quadratically with the number of refinement levels for 
our application. This is because, as described in the introduction, basis functions are kept 
in the basis whenever they contain an atomic nucleus within their support. All functions 
in this subset of significant size of the basis functions associated with each atomic center 
therefore touch one another, and the multiscale matrices contain dense blocks connecting 
all of these elements of a given center with one another. Because the number of functions 
associated with a given center grows linearly with the number of refinement scales /c, the 
number of operations required in the naive approach of multiplying directly by these dense 
matrices scales quadratically with the number of functions in the basis. For reference, a typ- 
ical number of refinement levels in electronic structure calculations of the lighter elements 
would be = 5, as employed in the carbon atom and the nitrogen molecule [0]. 

A comparison with the previous implementation in Fortran90 on the same processor 
(Superscalar SPARC Version 9, UltraSPARC) demonstrates the speed improvements and 
scaling which can be achieved with the new approach. The "time" axis is the CPU time 
taken by one application of the operator. The "k" axis represents the number of levels 
of refinement made in the basis and is proportional to the number of points in S. 

Figure 9 compares the runtimes of in three dimensions on a 1-level touching good 
basis with 3rd order interpolets consisting of concentric cubes of size 15'^ centered about one 
atomic nucleus, as would be appropriate for the calculation of the electronic structure of a 
single atom. Although there is initially a significant 0{N) contribution, as a function of the 
number of refinement levels k, the times for the naive approach show the constant increments 
in slope characteristic of a quadratic function. The new approach compares very favorably 
and is about thirty times faster for typical values of k. (Note the difference in vertical scale 
between the two figures.) 

Although the comparison in Figure 9 is quite favorable for the new algorithm, one must 
bear in mind that given the typical decay in the interpolet expansion coefficients about an 
atom0, |1|], the functions which are appropriate to maintain in the expansions tend to have 
the 2-level touching property, not the 1-level touching property. Figure 10 compares the 
runtimes of in three dimensions on a 2-level touching good basis of concentric cubes of 
size 9^, where the speed up just as dramatic as before, now by approximately a factor of 40. 

Figure 11 compares the runtimes of in three dimensions on a 2-level touching good 
basis of two refinement centers, with refinements now consisting of cubes of size 9^ (similar 
to figure 5). This situation arises in the the calculation of the electronic ground state of the 
nitrogen molecule, N2. Note that with the introduction now of two atomic centers the times 
are again consistent with the scalings described above: The runs times only double in the 
multilevel algorithm but quadruple in the naive algorithm. 

Having considered the efficiency of the algorithms, we next turn to the use of these 
algorithms in the solution of Poisson's equation to determine electrostatic fields, which was 
the rate limiting step in the calculations carried out in p[ and |]l|. From those calculations. 
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we were aware that the combination of conjugate gradients with the simple preconditioning 
consisting of just applying the inverse of the diagonal elements of the Laplacian matrix leads 
to an algorithm requiring very few iterations. It was the application of the operator within 
each conjugate gradient iteration which limited the efficiency of those earlier calculations. 

Figures 12-14 illustrate results for varying levels of refinement in the two different systems. 
The first system consists of refinements centered about a single atomic center within a 
cubic super cell of side 15 Bohr radii with periodic boundary conditions. (One Bohr radius 
is approximately 0.529 Angstroms.) The second system contains two refinement centers 
separated at a distance of 2 Bohr radii, approximately the inter-nuclear separation in the 
nitrogen molecule. This latter system resides within a rectangular supercell of dimensions 
(15 Bohr)^x(17 Bohr). In both cases, the spacing of the grid at the coarsest scale is 1 Bohr, 
and the finest spacing is Bohr. At k = 22, the greatest refinement considered in our 
numerical experiments, the finest grid spacing is approximately 0.24 x 10~^A. A full grid 
at this resolution would contain 2.8 x 10^^ points. Our truncated basis contains only about 
60,000 functions in this case. 

Figure 12 compares, as a function of the number of refinement levels k, the condition 
number of the Laplacian represented in a truncated interpolet basis (the "stiffness matrix" 
for the basis) and in an untruncated orthogonal basis at the corresponding resolution. The 
figure also shows the effect on the condition number of the interpolet stiffness matrix of the 
simple diagonal preconditioner described above. The condition numbers for the truncated 
interpolet bases were determined numerically using the operators implemented as described 
above. The curves indicate results for the system with a single atomic center, and the 
symbols indicate results for the two atom system. Comparing the results for the one and 
two atom cases suggests that apart from some transient behavior for small fc, the condition 
number is not sensitive to the number of atoms and depends primarily on the number of 
refinement levels k. 

Although finite basis representations of the Laplacian constructed from orthogonal func- 
tions at a given resolution should all have similar condition numbers, the fact that the 
interpolet basis is not orthogonal allows the condition numbers of multiscale interpolet oper- 
ators to be quite different that their single-scale counter parts. Compared to an orthogonal 
basis, the condition number in the interpolet representation is already over two orders of 
magnitude superior at the typical A; = 5 levels of refinement. This comparison continues 
to improve with increasing scale. The orthogonal condition number scales inversely as the 
square of the spacing on the grid of maximum resolution whereas the interpolet condition 
number scales inversely with approximately the 5/4 power of the resolution, as determined 
from the slope in the figure. The interpolet basis itself therefore provides an intrinsic form of 
preconditioning. Figure 12 shows that our simple explicit, diagonal preconditioner improves 
the scaling of the condition number, which now scales merely as the inverse of the resolution. 
(Note the lower slope of the lower curve.) At A; = 5 levels of refinement the improvement is 
only a factor of three but becomes more significant as the number of refinements increases. 
We extrapolate, based upon the observed scaling behavior, that the improvement is by a fac- 
tor of sixty aXk — 22 levels of refinement, the greatest refinement considered in the examples 
below. 

Figure 13 shows the convergence of the preconditioned conjugate gradient algorithm 
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in the solution of Poisson's equation. As a simple example, we solve for the electrostatic 
potential which arises from the two nuclei in a nitrogen molecule. In this calculation wc use 
k = 8 levels of refinement, a somewhat higher level of resolution than would be employed in 
calculations of the N2 molecule. For this calculation, the charge density of each nucleus is 
modeled as a three dimensional Gaussian of root mean square width a along each direction 
equal to the spacing on the finest scale. After an initial phase of about twenty iterations, the 
convergence becomes nearly perfectly exponential. This procedure reduces the magnitude of 
the residual vector by ten orders of magnitude in one hundred iterations. This is very good 
performance for a system consisting of 14,000 degrees of freedom with a Laplacian operator 
with a nominal single-scale condition number of about 65,000 at this level of resolution. 
The slope of this exponential portion of the convergence curve corresponds to a reduction 
in error at each iteration by 25%. One would obtain the same error reduction in a simple 
weighted iterative Jacobi algorithm (with the inverse of the maximum eigenvalue as the 
weight) applied to an operator with condition number c ~ 4. The quantity c, the inverse 
of the fractional improvement in the magnitude of the residual, we define as the effective 
condition number for the algorithm. 

Figure 14 shows this effective condition number c for the conjugate gradient algorithm 
with simple diagonal preconditioning as a function of the number of refinement levels k for 
the solution of Poisson's equation for the nuclei in the nitrogen molecule. In all cases the 
extent of the nuclei a is again set to the spacing of the finest grid. We note that after 
about six refinements, the effective condition number is essentially constant. The example 
from Figure 13 is therefore representative of the typical rate of convergence attained. These 
results indicate that, regardless of the amount of refinement, a constant number of iterations 
will suffice to produce a result of a given accuracy, even as the nominal condition number 
for an orthogonal stiffness at the corresponding resolution approaches 1.8 x 10^^ at A; = 22. 
Because the computational work involved in each iteration is linear in the number of points 
in the basis, this approach appears to produce the solution to Poisson's equation in 0(N) 
time for these multiresolution bases. 
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A Notation 



Ck 2^ ■ Z" for A; > 0, for k < 0. 

Dk Ck-i - Ck for A; > 0, for A; < 0. 

9k min(A;, m) where m is the largest integer such that 2™ divides all the components of y. 

Tk{S) the space of functions over S <Z Z'^. 

Xki(t>,S) the space of linear combinations of 0( ) for y e 5" C Z". 

it the mapping from S Ik{<P, S) which takes y (pi ), and linearly extended to a map fr( 

V the zero-lift representative oiv ^ , i.e. v e Tk-, such that {;(?/) = v{y),y e and v{y) — 0, ? 

y* the dual space of the vector space V. 

Jk^M the map, J^fc, Tk-^, given by Jk^^ = ° ''ki- 

Jk short for Jo k- 

J_fe short for Jk^. 

{f\0\g) the matrix element J f{x)Og{x)(r'x. 
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Figures: (available from [http:/ /laisla.mit.edu/mucliomas/Papers/nonstand- 



hgs.psl or in the file "figs.ps" in this posting) 

Figure 1: (a) shows 4>{x/2) = ^4>{x — 3) + ^0(a; — 1) + 0(a;) + + 1) + ^4>{x + 3). 

(b) shows the functions — 3), jQ<t>{x — 1), + 1), + 3). 

Figure 2: (a) the smooth function f{x) = (| + 3)^e~*^t)*. (b) the approximation to f{x) in 
Xq{Z). (c) component of the approximation in Xi(Ci). (d) component of the approximation 
in Jipi). 

Figure 3: Approximation of the function f(x) = e~'^' with interpolets centered at the 
tick marks: (a) a graph of e~'^'. (b) the cardinal approximation constructed from interpolets 
centered at the knots of a uniform grid, (c) and (d) cardinal approximations from local 
refinements of the uniform grid. The additional levels of refinement near the singularity 
increase the accuracy. 

Figure 4: (left) a graph of the relative error for k levels of local dyadic refinement, 
(right) a graph of the C°° relative error for k levels of local dyadic refinement. 

Figure 5: An example of a truncation one might wish to use for a diatomic molecule (two 
atomic cores). In black are shown those points in S G Z"- whose residual values might be 
significant. 

Figure 6: A visual summary of the 1-level touching condition. Solid plots represent 
functions centered at points in 5*. Dotted plots represent functions centered at points in 
— 5". Tick marks delimit the overlap of the two functions. 

Figure 7: Our example of the 1-level touching good basis in one dimension. Note that 
the two functions plotted do not touch. 

Figure 8: A generic example of a truncation which meets our definitions of good and 
1-level touching. 

Figure 9: The previously used implementation is on the left, and the implementation 
employing a 1-level touching algorithm is on the right. (Note the difference in scale on the 
vertical axes.) 

Figure 10: The previously used implementation is on the left, and the implementation 
employing a 2-level touching algorithm is on the right. (Note the difference in scale on the 
vertical axes.) 

Figure 11: The previously used implementation is on the left, and the implementation 
employing a multilevel algorithm is on the right. (Note the difference in scale on the vertical 
axes.) 

Figure 12: The condition number of the Laplacian operator represented in truncated 
multiresolution interpolet basis as a function of the number of refinement levels k with 
and without simple diagonal preconditioning and compared with the condition number in 
an orthogonal basis with the same resolution. Lines indicate results for bases with a single 
atomic center of refinement and points represent results for two atomic centers corresponding 
to the nitrogen molecule. 

Figure 13: Convergence of the solution to Poisson's equation for the nuclear potential in 
a nitrogen molecule in an interpolet basis with k = 8 levels of refinement about each nucleus. 

Figure 14: Effective condition number of Poisson's equation for the nuclear potential in 
a nitrogen molecule with simple diagonal preconditioning as a function of k the number of 
levels of refinement about each nucleus. 
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